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1 Introduction 

Studying cellular response to radiation damage is important from the perspectives of 
both radiotherapy and the mitigation of radiation damage. In the case of the former, 
the goal is the induction of cell death, whereas in the latter case a cellular response 
leading to organismal survival is required. 

The fate of the organism after radiation damage is linked in a complex manner to 
cellular fate. In the short term, an extensive cell death due to radiation damage will 
lead to the death of the organism. In the long term (relevant only if the organism 
survives the short term effects), it the survival of healthy and robust cells that dic- 
tate the organism's survival (for example, absence of long term cancerous cells). In 
this manuscript, we focus on the immediate survival of the organism after radiation 
damage. 

The mechanism of cell death after radiation damage is purported to be apoptosis 
via caspase activation that occurs several hours after radiation injury. Thus, it is 
natural to investigate cellular apoptotic machinery. Cellular response to radiation is 
very complex and involves, possibly, proteins that respond to DNA damage and the 
formation of free radicals that modify cellular biochemistry. Although the exposure to 
radiation may be momentary, the effect of the exposure on cellular biochemistry may 
be long lived. Further, several proteins that are expressed transiently after radiation 
damage, may trigger downstream response such as caspase activity that occurs long 
after their expression. 

In view of the complex cellular response, protein-protein interactions that lead 
to that response, and possible importance of temporal evolution of the system, it is 
important to study systemically the time-dependent cellular response to stress. 
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There have been various efforts to model cellular response to radiation, as well as 
to model apoptosis. Efforts to model cellular response to radiation were largely pre- 
cipitated by the observation, at the level of a single cell, of oscillatory p53 response to 
radiation damage. 1 In particular, the oscillatory response of p53 to radiation damage 
has been modeled mostly using deterministic simulations with kinetic rate laws,- - - 
and, to a significantly lesser extent, via stochastic simulations.- Cellular apoptosis 
has been modeled, independent of p53 radiation response, via deterministic simula- 
tions by several research groups,- - - including efforts to include cell-cell variability in 
a probabilistic manner.—^ 

Recently, there have been a few efforts to link the p53 response to cell sur- 
vival/death using deterministic simulations,— _ — or with inclusion of limited stochas- 
ticity. 15 Due to increasingly larger pool of experimental observation on the proteins 
involved, and their cellular compartmentalization, in the apoptotic network, model- 
ing of transient system behavior remains an active area of research and serves as a 
platform for evaluating pharmacological strategies. 

Spatio-temporal modeling using stochastic simulation methods allows for detailed 
examination of the kinetics of biochemical reaction networks within a cell. From 
such simulations, one can address several different issues: spatial localization of reac- 
tions (e.g., on the mitochondrial membrane), as well as studying the time evolution 
of molecular concentrations in response to perturbations to the system. These per- 
turbations can be radiation stress on the cell, or the response to a given drug dose. 
Further, stochasticity allows for the inclusion of cell variability in a natural way, as 
well as for studying conditions where the number of certain types of molecules become 
very low and their presence/absence dictates system behavior. 

In this study, we perform stochastic simulations to study the time evolution of pro- 
teins, associated with apoptotic pathways, in response to radiation damage. Further, 
we check the efficacy of various treatment strategies, including polypharmacological 
effects. In particular, we use different inhibitors targeting different proteins in the 
reaction network. The manuscript is organized as follows. We first present the bio- 
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chemical reaction network, gathered from experimental observations, associated with 
radiation induced apoptosis. Then, we describe briefly the stochastic simulations 
that we perform to study cell response to radiation damage, followed by results and 
discussion. 

2 Biochemical reaction network 

We present here the chain of biochemical events that we use to model apoptotic 
mechanisms that are activated after radiation damage. Instead of adopting a com- 
prehensive model that includes all the cellular processes that are activated, or affected, 
by radiation damage, we focus on a number of key proteins/genes as a first approxi- 
mation. An overview of the key interactions is given in Figure [TJ Tabled] gives details 
about the proteins involved in the system. Several other proteins besides those in 
Figure [1] are included in the model and described in Table [TJ since an interaction 
map, such as that in Figure (TJ does not contain details of the biochemical events that 
lead to the individual interactions. 

Before we describe the radiation response of the network, we present the generic 
underlying network. 

2.1 p53 module 

p53 is an integral component of cell death due to genotoxic stress (including ra- 
diation damage), and it activates apoptosis via both transcription-dependent and 
transcription-independent pathways.— _ — The exact effect of p53 in its transcription- 
independent role in apoptosis is still in-debate (see a recent review by Lindenboim et 
al.—), we attempt to incorporate some of the key experimental observations. 

The response of p53 to radiation exposure has been studied in great detail at 
the individual cell level. The response of p53 to radiation determined experimentally 
forms the starting point of our modeling effort: we first focus on the p53 "module" 
(as highlighted in Figure [1] with the dashed rectangle) and calibrate our model to 
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reproduce the known p53 response. 

p53 exists in several forms in a cell based on its post-translational modifications 
(e.g., phosphorylation, ubiquitination at different sites). In this manuscript, we utilize 
a model that captures the known biological features: such as p53 oscillatory behavior 
upon radiation damage,- 1 ^^ and a stress-induced increase in nuclear p53^i and its 
transcriptional activity.— 1 ^ 

p53 is translated in the cytoplasm (p53C) and this formation is represented by 

$ ^ P 53C (1) 

where $ represents a null state, and the above reaction represents the formation of 
p53C under normal operating conditions given by the rate K\. 

p53C can translocate to the nucleus, especially when a cell undergoes stress.— 
This translocation is modeled by 

P 53C P 53n (2) 

where p53n is the nuclear p53. In the nucleus, p53n transcriptionally activates Mdm2 
mRNA and this transcription activation is modeled as cooperative:- p53n tetramers 
on the DNA leads to this transcription activation. The oligomerization of p53n is 
depicted by 

4p53n p53no (3) 

fc 3 

where the reaction order of the forward reaction is 4. The oligomerized p53 leads to 
a direct transcription activation of Mdm2: 

p53no — ^» p53no + Mdm2m (4) 

The Mdm2 messenger RNA translocates to the cytoplasm where it is translated. 
The cytoplasmic Mdm2 (Mdm2C) can be imported into the nucleus. This series of 
events is described by the following reactions: 

Mdm2m ^ Mdm2mC (5) 
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Mdm2mC ^ Mdm2C (6) 

Mdm2C Mdm2 (7) 

The monoubiquitinated, nuclear p53 (p53NU) is formed via the interaction of 
p53n with Mdm2 (or similar ligases),— and we represent this by 

k + 

p53n + Mmd2 P 53.Mdm2 ^ p53NU + Mdm2 (8) 

k s 

A similar reaction is modeled to occur for the cytoplasmic p53 and Mdm2: 

k+ 

p53C + Mmd2C ^ p53C.Mdm2 p53CU + Mdm2C (9) 

kg 

Further, the ubiquitinated nuclear p53 can translocate to the cytoplasm: 

P 53NU p53CU (10) 

The ubiquitinated cytoplasmic p53 can translocate to the mitochondria^ as shown 
below: 

p53CU ^> P 53m (11) 

Upon severe genotoxic stress, p53n become transcriptionally active for pro-apoptotic 
proteins^ such as PUMA and Bax, and such a transcriptionally active form is rep- 
resented by p53T: 

p53n ^ p53T (12) 

Although the major role of p53T is transcriptional activation of pro-apoptotic pro- 
teins, we let it retain the capacity for oligomerizing and activating Mdm2 (this, 
though, is a minor effect), 

4p53T v=f p53To (13) 

^13 

p53To ^ p53To + Mdm2 (14) 

The increase in p53 upon radiation damage is not expected to be a step function. 
Similarly, return of p53 to baseline levels after removal of radiation source is not 
expected to be instantaneous. Indeed, oscillatory levels of p53 persist for an extended 



period of time after radiation damage.- Further, p53 response to radiation is mediated 
by upstream proteins such as ATM. We do not model such upstream agents in specific 
detail, instead we represent the upstream events by the following set of reactions. 



$ ^ A (15) 

^15 

A^>A + B (16) 

B ^ $ (17) 

B^%B + p53C (18) 



Thus, A activates B, which in turn activates p53C. The effect of these upstream 
reactions is a smooth increase in p53C concentration after radiation damage. 

In addition to the above equations, we model the independent degradation of 
p53C, p53n, p53CU, p53NU, p53T, and p53m to highlight the fact that external 
agents (i.e., those not modeled explicitly) have a role in modulating protein concen- 
trations. For example, several E3-ligases (other than Mdm2) ubiquitinate and label 
p53 for proteasomal degradation downstream. These reactions are given below: 



p53C ^ $ (19) 

P 53n ^ $ (20) 

p53CU $ (21) 

p53NU ^ $ (22) 

p53T ^ $ (23) 

p53m ^4 $ (24) 



2.2 p53 modulated apoptotic pathways 

Now, we discuss the mitochondrial apoptotic pathways downstream of p53, and the 
role of p53 module to the activation of mitochondrial apoptosis (see Figured]). 
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p53T leads to transcription activation of proapoptotic PUMA and Bax, as modeled 

by 

p53T ^ p 53T + PUMA (25) 

and, 

p53T ^% p53T + Bax (26) 

The roles of PUMA and Bax are discussed below. 

p53m binds to the anti-apoptotic protein Bcl-2 to inhibit the anti-apoptotic ac- 
tion of Bcl-2 on the mitochondria.— The biochemical equation associated with this 
reaction is given by 

27 

p53m + Bcl-2 p53m.Bcl-2 (27) 

_>7 

with binding affinity given by k^/k^. PUMA also binds to the antiapoptotic Bcl-2 
member on the mitochondrial membrane. 27 

k + 

PUMA + Bcl-2 ^ PUMA.Bcl-2 (28) 

K 28 

Further, PUMA binds stronger to Bcl-2 than does p53, and can displace p53 from 
its complex with Bcl-2— - as can be seen by a combination of reactions |27] and EHJ 
thus, freeing up p53 for further proapoptotic activity on the mitochondria, such as 
by activating Bax/Bak as shown below. 

Pro-apoptotic Bax translocates back and forth from the cytoplasm to the mito- 
chondria 28 to establish an equilibrium between the cytosolic and mitochondrial Bax 
concentrations: 

k+ 

29 

Bax Bax mito (29) 

At the mitochondria, Bax interacts with Bcl-2 to form a complex— that prevents Bax 
oligomerization: 

k+ 

3Q 

Bax mito + Bcl-2 Bax mito . Bcl-2 (30) 

k~ 

Thus, the available Bcl-2 (that also forms complexes with p53 and PUMA as shown 
above) regulates the amount of Bax available for oligomerization. 
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Chipuk et al. 27 reported that p53m binds stronger to Bcl-2 than does Bax, and 
can displace Bax from its complex with Bcl-2. This reaction can be modeled via the 
following reaction: 

fc + 

p53m + Bax mito . Bcl-2 p53m. Bcl-2 + Bax mito (31) 

^31 

However, eq ED is a straightforward linear combination of eqs [27] and [301 an d, thus, 
the obvious requirement of consistent equilibrium constants (k^/k^ = k^k^/ k^ 7 k^ ) 
must be observed (or, eq [31] should not be considered explicitly). 

Active form of Bax on the mitochondrial membrane is needed for oligomeriza- 
tion, a process that requires further interaction of Bax with activators such as tBid. 
tBid localization to mitochondria^ 3 ^ is an important event for Bax activation and is 
modeled by the following reaction: 

k + 

tBid ^ tBid mito (32) 

^32 

Subsequently, tBid activates mitochondrial Bax into an active (for oligomerization) 
form: 

tBid m i t o + Bax mit o tBid m i t o + Bax* (33) 

It has also been suggested that tBid activates Bax in the cytosol which then localizes 
to the mitochondria. 31 However, the order of consecutive reactions is not expected to 
have a significant effect further downstream (such as caspase activation). A similar 
Bax (or Bak) activating effect can be achieved by p53m:— 

p53m + Baxmito — ^> p53m + Bax* (34) 

and by PUMA, 32 

PUMA + Baxmito ^ PUMA + Bax* (35) 

The activated Bax then oligomerizes on the mitochondrial membrane: 

k + 

Bax* + Bax* ^ Bax 2 (36) 
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Further, higher order oligomerization is modeled following Albeck et al.— 

37 

Bax 2 + Bax 2 Bax 4 (37) 

^37 

Mitochondrial outer membrane permeabilization (MOMP) induced by Bax oligomer- 
ization results in the release of cytochrome- c from the mitochondria into the cyto- 
plasm. We model the release of cytochrome- c via the following process:- 

Bax 4 + cytc mito — ^ Bax 4 + cytc (38) 

where the above reaction serves to model the fact that the presence of a Bax pore 
is essential for release of cytochrome- c into the cytoplasm. Similarly, MOMP also 
releases Smac/Diablo into the cytoplasm: these proteins bind to anti-apoptotic IAP's 
and, thus, facilitate apoptosis. 33 

Bax 4 + Smac mito — Bax 4 + Smac (39) 

Cytoplasmic cytochrome- c interacts with Apaf in an ATP-dependent manner, and 
the heterodimer forms a heptameric complex called the "apoptosome".— Following 
Bagci et al.—, we model these events via the following two reactions: 

k + 

cytC + Apaf v=^ cytC.Apaf (40) 

7cytC.Apaf apop (41) 

Further, we use cooperativity in the formation of the apoptosome: the order of the 
forward reaction |4T] is 4. 

The apoptosome then incorporates caspase-9, that in its autocatalyzed, activated 
form cleaves procaspase-3 (C3) to form active caspase-3 (C3*). Following Albeck et 
al.,— we use the following reactions for these processes: 

fc+ 

apop + C9 ^ apop.C9 (42) 



apop.C9 + C3 ^ apop.C3.C9 ^ apop.c9 + C3* (43) 



As mentioned above, tBid is an important factor for activation of the caspase 
cascade. Caspase-8, activated due to external death signals, is an important molecule 
responsible for Bid cleavage." Additionally, active caspase-3 also results in Bid cleav- 
age, resulting in a positive feedback loop activating apoptosis. 36 We model the asso- 
ciated reactions following Bagci et. al. 7 : 

C3* + Bid ^ C3 * .Bid ^4 C3* + tBid (44) 

^'44 

k + 

C8 + Bid ^ C8.Bid C8 + tBid (45) 

The concentration of active caspase-8 depends on external death factors. 

XIAP inhibits the activity of apoptosome^2& and promotes proteasomal degra- 
dation of caspase-3— preventing cell death. On the other hand, Smac inhibits the 
activity of XIAP.— These processes are modeled as,- 

k + 

apop.C9 + XIAP apop.C9.XIAP (46) 

k 1 r, 



C3* + XIAP ^ C3*.XIAP ^ C3* Vh + XIAP (47) 
C9 + XIAP ^ C9.XIAP ^ C9 Ub + XIAP (48) 



K 48, 

h 

h 4<\ 



Smac + XIAP v=f Smac.XIAP (49) 

In addition to above reactions, there are formation reactions for several species - 
Bcl2, Bax, Apaf, C9, C3, XIAP, Smacm, cytcm, C8, and Bid. Further, these species, 
along with C3*, C3U, C9U, Bax*, and p53 forms mentioned above also undergo 
degradation reactions to help establish steady state. 

3 Drug interactions 

In this section, we integrate the possible drugs to regulate apoptosis in normal cells 
into the biochemical network discussed above. In particular, we consider four specific 
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proteins targets - Mdm2, PUMA, Bid, and C3* (Mdm2-I, PUMA-I, Bid-I, and C3-I, 
respectively). The actions of the inhibitors are modeled by the following reactions: 

Mdm2.Mdm2 - I, (50) 

PUMA.PUMA - 1, (51) 
Bid.Bid - I, (52) 

C3*.C3 - I. (53) 

Known inhibitors of p53/Mdm2 interaction includes nutlins, BEB55, BEB59, 
and BEB69. Mustata et al.— have identified several PUMA inhibitors using phar- 
macophore modeling, and several of these compounds have been shown to inhibit 
PUMA-induced apoptosis in vitro. 

4 Stochastic simulation method 

The number of proteins in a typical cell (about picoliter) are significantly less than 
Avogadro's number: a nanomolar concentration of a protein corresponds to ~ 600 
molecules of that protein in the cell. Accordingly, for smaller cells and/or lower con- 
centrations, deterministic rate equations are not applicable. For this reason, we use 
the well-known Gillespie algorithm for stochastic simulations of a system of chemical 
reactions. 

Details of the Gillespie Algorithm are given in the original papers.— ^2 Here we 
briefly describe some of the ideas behind the algorithm. 

For illustration, consider the following schematic reaction, 

A + B -> AB, (54) 

and denote the stochastic rate constant by c. This stochastic rate is directly related 
to the macroscopic kinetic rate, k, by a simple relation (which, for the reaction in eq[T] 
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and 



K r,o 



Mdm2 + Mdm2 - I 

k l 

PUMA + PUMA - I ^ 

k 



51 
k+ 



Bid + Bid -I 

k 



■12 



C3* + C3 - I 



is c = kV, where V is the reaction volume). At any instant of time, the propensity 
for the reaction in eq [52] is 

a a = N A N B c (55) 

where and are the number of molecules of A and B in the reaction volume, 
and the subscript a denotes that the reaction index is a in the system of M chemical 
reactions (1 < a < M). 

In the Gillespie algorithm, one such reaction is chosen to be proportional to its 
propensity (we use the so-called Direct Reaction version), and the time is advanced 
based on the overall reaction propensity at that time.— & The underlying assumptions 
are that the system is in thermal equilibrium (the number of reactive collisions are 
much less than the number of non- reactive collisions), and that it is well-mixed. 

4.1 Heterogeneity 

Real cellular systems are heterogeneous - e.g., mitochondria and cytoplasm offer vastly 
different environments (and, there may be heterogeneities within these compartments 
themselves). Thus, at a first glance, the well-mixed assumption of the Gillespie algo- 
rithm seems to preclude an application to real cellular systems. However, Bernstein- 
developed a protocol to extend the Gillespie algorithm to heterogeneous environments. 

In that protocol,— the cell (or, any simulation volume) is subdivided into elements 
and a reaction is reproduced and treated as an individual reaction in each subvolume. 
Accordingly, each reactant is treated as a different molecule in each subvolume {e.g., 
Aj 7^ Aj, where i and j are different subvolumes), and the reaction propensities 
are altered accordingly. Further, each subvolume is well-mixed - thus allowing for 
formulating the system using the Gillespie algorithm. In context of this work, the 
nucleus, the cytoplasm, and the mitochondrial membrane are well-mixed, distinct 
subvolumes. 

This protocol results in an increase in the number of reactions both due to the fact 
that each reaction is treated independently in different subvolumes, as well as addi- 
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tional "reactions" (due to diffusion) relating to the conversion of A* to Aj, as described 
in detail by Bernstein.— Based on the identity of a compartment/subvolume, a type 
of protein may not be present in that compartment (for example, the apoptosome 
does not permeate through the mitochondrial membrane). 

5 Model parameters 

We obtain several parameters from experimental results, and obtain others in the 
range of previous computational studies.-^ Here we discuss some known experimental 
rate/equilibrium constants, and the initial concentrations of the species involved in 
the biochemical reaction network discussed above. 

The dissociation constant, k^j/k^, for p53 and pro-apoptotic Bcl-2 family mem- 
bers have been given variously as 160 nM— and 535 nM— - both obtained via surface 
plasmon resonance. We choose an intermediate value of the dissociation constant 
(333 nM). PUMA binds significantly more strongly, and Chipuk et al— report the 
dissociation constant k^g/k^g = 10 nM, a value that we use. 

Edlich et al. 28 report the on and off rate of Bax from cytosol and mitochondria 
as approximately the same, and establish a first order reaction (as represented by 
reaction E]) with fc+ = 4.9 x 10~ 3 s" 1 and = 4.7 x 10~ 3 s -1 . 

For reaction (3TJ it is reported that p53 displaces Bax from its complex with Bcl-2 
at equimolar concentration, whereas Bax displaces p53 from its complex with Bcl-2 
at a 50 times higher concentration. 18 This implies that fc^/fc^ = 50. Since eqs [271 
1301 and EH are not independent, the above stated values imply that k^ /k^ ~ 17 /jM. 
This value of dissociation constant is significantly higher than values reported in the 
literature (0.1 /xM— and 20 nM— ). We use kJ /k^ Q nM, except where noted explicitly 
- and, as shown below, the results obtained are robust to a wide range of values for 
this dissociation constant. 

Half lifes of a few proteins have also been reported in vivo. However, the direct use 
of such an information (for example, to compute the decay rates of the type A— >• $) 



13 



is often not appropriate because the decay rate of protein A in vivo is convoluted 
with, for example, the decay of protein B that upregulates A. 

The parameters for the p53 module that differ for these three levels of radiation are 
given in Table [21 and the parameters that stay the same irrespective of the radiation 
level are given in Table |3j These parameters for the p53 module are chosen to reflect 
p53 oscillations observed experimentally. Ku and K12 do not affect the dose-dependent 
oscillatory behavior of p53, however, they are important for coupling of the p53 
module to the apoptotic network, and the corresponding values are given in the text 
and figures below. 

5.1 Initial concentrations 

When a normally unstressed cell is exposed to radiation, it undergoes a sequence of 
biochemical events that may, depending upon the radiation dose, lead to apoptosis. 
Thus, the appropriate initial conditions of cellular protein levels for such an apoptotic 
study are the steady state conditions in an unstressed cell. 

Irrespective of what initial number of molecules we start with, we first allow the 
system to establish a normal steady state (without radiation-induced apoptosis) by 
selecting N-IR kinetic parameters from Table [2] (in addition to radiation-independent 
kinetic equations/parameters). Subsequently, we perturb the system with radiation 
damage (modeled here as genotoxic effect), and allow the system to evolve. With 
enough stress, the system will eventually display apoptosis, with a high concentra- 
tion of pro-apoptotic proteins (such as active caspase-3). In essence, the initial con- 
centrations that are used to study radiation induced damage and subsequent drug 
treatments are the steady state conditions without any radiation/drug. 

For reference, the typical concentration a typical protein present in the cell under 
normal conditions is 10-100 nm.- In a picoliter cell, the number of molecules of a 
typical protein is, thus, in the range 6000-60000. On the other hand, proteins such as 
active caspase-3 (C3*) are expected to be low under normal homeostatic conditions. 
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5.2 Modeling radiation damage 

Mitochondrial apoptosis is regulated by transcription-dependent and independent 
role of p53. Thus, we focus on the p53 module (Figure [T]) to model the effects of 
radiation damage. The downstream action of the apoptotic pathways (for example, 
represented by Bax oligomerization, cleavage of caspase-3) occur due to the response 
of the p53 module to radiation damage. Specifically, we model the response of p53 to 
three different radiation levels: no radiation (N-IR), low radiation (L-IR), and high 
radiation dose (H-IR). 

The procedure for the simulation is as follows. (1) Establish steady state for an 
unstressed cell (N-IR). (2) induce radiation damage (L-IR or H-IR), as modeled by 
use of the appropriate kinetic parameters from Table [2] for a certain time AT. (3) 
Return to unstressed parameters after AT. 

We note here that AT does not refer to the time duration of radiation exposure, 
rather it refers to the amount of time the radiation modifies the kinetic parameters. 
AT only indirectly refers to the exposure time: longer exposure to a given radiation 
dose will affect cellular biochemical kinetics for a longer duration. 

6 Results 

6.1 Radiation induced p53 oscillations 

Single-cell experiments to monitor cellular p53 levels after radiation damage showed 
that individual cells show dramatic oscillations in p53 levels. In particular, these 
oscillations are sustained and show dose-dependent oscillation frequencies.- These 
experimentally observed results form the basis of simulating the p53 module using 
stochastic simulations: we reproduce the radiation dependent behavior of p53 re- 
sponse using stochastic simulations. 

Figure [2] shows the oscillations in p53n (the major p53 component upon radiation 
damage) as a function of time for two different radiation doses and exposures. The 
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chosen model for radiation damage reproduces the known p53 oscillations, that are 
dependent upon the radiation dose. At low levels, this frequency is ~ 11 hr, whereas 
it decreases to w 6 hr upon increasing the radiation dose. 

The differences in L-IR and H-IR parameters for the p53 module suggests one 
possible mechanisms for the observed radiation-dose dependent frequencies. Impor- 
tant factors in the model leading to this difference are the enhanced translocation of 
Mdm2 protein from the cytoplasm to the nucleus, and the enhanced translocation of 
Mdm2 mRNA from the nucleus to the cytoplasm. This was deliberately chosen to 
represent the time delay in Mdm2 response that is frequently utilized in deterministic 
models.- 1 ^ 1 ^ 

Thus, this study shows that such time delay can be modeled mechanistically in 
stochastic simulations (and allow for the modulation of oscillation frequency) by the 
realization that proteins/genes act in specific cellular contexts. Further, the increase 
in the translocation rate of the Mdm protein (along with increased translocation of 
p53) to the nucleus upon DNA damage has been suggested.— 

Further, as Figure [2] suggests, at low exposure time (modeled by AT = 2 x 10 4 s), 
there is a single peak in the p53 profile for both radiation levels. On the other hand, 
sustained p53 oscillations are obtained for high exposure times (AT = 2 x 10 5 s). Even 
for cells exposed to radiation for the same amount of time, it is possible that cell-cell 
variability can lead to different amounts of observed oscillatory peaks (AT is, after 
all, not the radiation exposure time, but is the time for which the kinetic parameters 
are altered from unstressed levels). 

6.2 Role of p53 oscillations in apoptosis 

At higher radiation dosage, the oscillations in p53 occur at a higher frequency. How- 
ever, it is unclear if cellular apoptotic response is mainly determined by the oscillation 
frequency: is increased apoptosis at higher dose of radiation a direct result of the 
higher p53 oscillation frequency? 

To address this question, we compare the caspase 3 activity obtained via the use 
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of H-IR and L-IR parameters in Table [2j along with the kinetic parameters associated 
with the rest of the apoptotic machinery. Further, we use the same rates for apoptosis- 
related transcriptional activation of p53 (k%2 = 1CT 5 s _1 ), and for p53 translocation 
to mitochondrial membrane (ft n = 10~ 5 s _1 ). The obtained caspase-3 concentrations 
for a cell are shown in Figure |3] for both H-IR and L-IR sets. 

Upon using a radiation-dose independent coupling of the p53 module with the 
apoptotic network, L-IR leads to an increase in caspase-3 activity upon radiation 
damage (radiation damage occurs at t = 0). On the other hand, H-IR leads to low 
caspase-3 concentration. Clearly, this result is at odds with experimental observation 
- increased radiation dose leads to an increase in apoptosis. Correspondingly, if K n 
and/or K12 are increased for H-IR, caspase-3 activity is increased (shown in Figure |3]). 

A crucial insight that emerges from Figure |3] is that the dose-dependent radiation 
damage and cell death depend upon the coupling of the regulatory p53 module to the 
apoptotic machinery (defined via an increase in the p53 pro-apoptotic transcriptional 
activity, K42, and an increase in p53 translocation to the mitochondria, K n ), and not 
merely on the p53-oscillation frequencies. We use the term "coupling" to describe 
how the p53-regulatory module leads to apoptosis via transcription-dependent (/C12) 
and transcription-independent («n) pathways. 

6.3 Downstream apoptosis cascade 

The response of the apoptotic network downstream of the p53 module sheds further 
light on the role of pro- and anti-apoptotic proteins in the cell. Here, we discuss 
(i) oscillations downstream of p53, (ii) role of Bax-Bcl-2 interactions, and (iii) the 
requirement for sufficient Bid cleavage for sustained apoptosis. 

Figure H] shows the concentrations of PUMA, Bax*, tBid, and caspase-3 after 
radiation damage with high dose (the qualitative response of low radiation dose is 
similar). 

PUMA mirrors the p53 oscillations, irrespective of the coupling strength. This 
phenomenon is due to direct activation of PUMA by p53T. Unlike oscillations in 
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PUMA level, oscillations in the concentration of activated Bax show a distinct de- 
pendence on the coupling strength: at low coupling strength (no apoptosis), Bax 
oscillations are sustained, but decay along with PUMA oscillations. On the other 
hand, a high coupling that results in apoptosis leads to a sustained increase in ac- 
tivated Bax at the expense of oscillations. Bax* increases initially due to the direct 
action of PUMA on Bax. This leads to the release of caspase-3 that cleaves tBid 
resulting in further Bax activation and a positive feedback loop. 

A sustained caspase-3 activity even after the decay of PUMA depends upon 
whether sufficient tBid is activated by the time PUMA decays. A sustained release 
of tBid and the positive feedback is very important for sustained caspase activity: 
for low coupling of p53 module and apoptotic network, insufficient PUMA activity 
results in an insufficient amount of caspase 3 release and tBid activation to maintain 
a sustained caspase activity. 

6.3.1 Strength of interaction between Bax and Bcl-2 

In Section [5], we discussed the strength of Baxm and Bcl-2 interactions ( ^30 / ^30 ) : dif- 
ferent values of this dissociation constant (either directly or indirectly) are suggested 
by different groups spanning several orders of magnitude.— Here, we discuss the 
effect of a range of this dissociation constant on the apoptosis model we use in this 
manuscript. 

Figure [5] shows caspase-3 and Baxm.Bcl-2 concentrations for two different values 
of the dissociation constants. A comparison of panels (a) and (c) shows that the 
downstream apoptotic activity is unaffected by dissociation constants that differ by 
two orders of magnitude (all other model parameters remain unchanged between these 
two panels): caspase-3 concentration is unaffected, although vastly differing amounts 
of Baxm. Bcl-2 complex are formed. Similarly, the different dissociation constants had 
an insignificant effect on systems showing a lack of apoptosis (panels (b) and (d)). 

The model is, thus, robust with respect to the strength of Baxm and Bcl-2 in- 
teraction. This suggests that cells can exhibit significantly differing strengths of this 
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interaction without affecting the cell fate given a genotoxic insult, and several alter- 
nate mechanisms can lead to the same cell fate. 

6.4 p53 transcription-independent apoptosis 

So far, we have discussed p53 transcription dependent apoptosis via transcriptional 
upregulation of PUMA and Bax upon radiation damage. On the other hand, the role 
of p53 in a transcriptionally independent manner is also of importance. In its tran- 
scription independent role, it has been suggested that p53 in the cytoplasm can either 
directly activate Bax, or can indirectly activate Bax by binding to anti-apoptotic Bcl- 
2 and, thus, preventing the latter's anti-apoptotic action. In this section, we explore 
this issue. 

Figure [6] shows the evolution of p53m (top panel) and caspase-3 (bottom panel) 
with time upon exposure to H-IR for different ratios of transcription dependent (/C12) 
and independent (acu) apoptosis pathway strengths. In the figure, the black line 
represents the case where the strength of transcription dependent pathways alone is 
not sufficient to cause apoptosis. Upon an increase in k u (black, red, and blue lines 
- in that order), p53m shows an increasingly upregulated behavior. Correspondingly, 
caspase-3 shows a sustained activity with an increase in p53 transcription independent 
pathway strength (for the same K12). 

Strikingly, the removal of direct Bax activation by mitochondrial p53 (green line) 
results in complete absence of apoptosis for the p53 transcription-independent apop- 
tosis, even for a high strength of transcription independent pathway. This suggests 
that the indirect role of p53m in apoptosis by binding Bcl-2 is not sufficient for sus- 
tained apoptosis by the p53-transcription independent mode - direct Bax activation 
by the mitochondrial p53 appears to be critical, too. 
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6.5 Role of inhibitors of specific proteins 

From a pharmacological viewpoint, it is important to identify protein targets that are 
good candidates for drug treatment to mitigate radiation damage. In this section, we 
focus on this issue with the aim of identifying potential targets suitable for treatments 
that are effective even if not administered immediately. 

Figure [7J illustrates the efficacy of treatments administered after two different de- 
lays: treatments administered immediately after radiation damage (15 minute delay, 
top panel), and treatments administered after a longer delay (12 hours, bottom panel). 
The radiation dose in Figure [7] leads to sustained caspase-3 activity in absence of drug 
treatment (bottom panel of Figure IU blue curve) . 

If treatment is available immediately after radiation damage, inhibition of PUMA 
is very effective in mitigating radiation damage. Similarly, inhibition of Bid also helps 
in mitigating radiation damage (although to a lesser extent than PUMA inhibition). 
In contrast, treatment via inhibition of Mdm2 and caspase-3 are not effective. On the 
other hand, if there is a longer delay before treatment is administered, PUMA inhibi- 
tion is an ineffective treatment and inhibition of Bid is the most effective treatment. 
The ineffectiveness of PUMA inhibition in this case results because PUMA activity 
till the administration of PUMA-I leads to sufficient activation of positive feedback 
loop involving caspase-3 and Bid: activation of Bax via tBid dominates the apoptotic 
response at this late stage. 

The main cause of apoptosis for the system in Figure [7J is the p53-transcription 
dependent pathway. Thus, Mdm2 inhibition, that decreases the ubiquitination of 
p53n and, subsequently, leads to an increase in p53T (and PUMA and Bax), increases 
the propensity of the system to undergo apoptosis. A comparison of the two panels 
of Figure [7J shows that administering Mdm2-I sooner after radiation damage leads to 
a quicker increase in caspase-3 activity. 
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7 Discussion 



7.1 Testable hypotheses 

Several key, testable features of cellular response to apoptosis emerges from the model 
and are listed as following, (i) oscillatory response of PUMA, at the level of a single 
cell, to radiation damage, (ii) elevated tBid activity even after PUMA decays when an 
increase in apoptosis occurs, and (iii) inhibition of Bid is more effective in mitigating 
radiation damage than PUMA inhibition when the treatment is administered after a 
substantial delay. 

The last point is especially relevant with respect to developing a treatment strategy 
for mitigating radiation damage. Even a very potent inhibitor of PUMA may not 
prevent apoptosis when it is administered after a substantial delay. Indeed, the 
inhibitor used in the current model acts in nM concentrations, and despite being 
effective when administered immediately after radiation damage, is ineffective after a 
longer delay because of enough tBid activity to maintain a sustained caspase activity. 

Even if potent inhibitors of Bid are not currently known, it is possible to test the 
last hypothesis above using the following two approaches. Firstly, a si-RNA knockout 
of PUMA performed several hours after radiation damage should be ineffective in 
mitigating apoptosis due to radiation damage. Secondly, knocking out Bid several 
hours after radiation damage should be more effective in mitigating apoptosis due to 
radiation damage. 
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Figure 1: A map of interactions in the apoptotic pathway. This map only shows a 
schematic of the interactions highlighting the main features of the model. Detailed 
biochemical reactions (that may involve proteins not specifically mentioned in this 
map) are given the text. 
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Figure 2: Sustained oscillations in p53 concentration as a result of the duration of 
alteration, via radiation, of cellular biochemical kinetics. Radiation damage at t — 
leads to radiation-dose dependent p53 oscillations. Single spikes in p53 concentrations 
(the left two panels) are observed when radiation induces a momentary change in 
cellular biochemical kinetics, and sustained oscillations (the two panels on the right), 
occur when the cellular biochemical kinetics are altered for a longer duration. Further, 
the frequencies of oscillations at the two radiation dosage (top right and bottom right 
panel) are similar to the experimentally observed values. 
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Figure 3: An increase in pro-apoptotic transcriptional activity of p53 associates with 
higher radiation dose to induce apoptosis. Caspase activities induced by the higher 
radiation dose for two different pro-apoptotic transcriptional activity of p53 (k^)- A 
lower dose of radiation leads to a substantial caspase activity at a smaller value of 
«42- K n — 10~ 5 for the three cases. 
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Figure 4: Role of tBid activation for sustained caspase-3 activity. A comparison of 
the concentrations of PUMA, caspase-3, and tBid obtained after irradiation with a 
high dose at t — with two different couplings (^12) - low value of K12 in the top 
panel does not lead to enough tBid activity to sustain apoptotic activity k u = 1CT 5 
for both panels. 
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Figure 5: Robustness of the model (for caspase-3 activity) with respect to the strength 
of Baxm and Bcl-2 interactions. Panels (a) and (c) show the results for k 12 = 2 x 
10~ 5 s _1 showing distinct apoptosis despite vastly different strengths of Baxm and 
Bcl-2 interaction. Panels (b) and (d) show an absence of apoptosis for different 
strengths of Baxm-Bcl-2 interactions (and both are for k 12 = 10~ 5 s _1 . 
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Figure 6: Important role of direct Bax activation by mitochondrial p53 in inducing 
apoptosis via the p53 transcription-independent pathway. The top panel shows the 
mitochondrial p53 concentrations for different sets of («ii,/ci2) values, and the asterix 
for the green curve indicates that direct Bax activation by the mitochondrial p53 in 
this case is abrogated (k 34 = 0). The bottom panel shows the corresponding C3* 
concentrations. 
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Figure 7: Effective mitigation of radiation damage via Bid inhibition at longer times. 
The top panel shows caspase-3 activity upon treatment via Mdm2, PUMA, caspase- 
3, and Bid inhibition after treatment by the respective inhibitors 15 minutes after 
radiation damage, and the bottom panel shows caspase-3 activity after treatment 
administered 12 hours after radiation damage. In the top panel, caspase-3 activity 
after treatment with PUMA inhibitors (PUMA-I) is indistinguishable from zero. In 
the bottom panel, resulting caspase-3 activity after treatment by PUMA-I (red line) 
is similar to the activity after treatment by procaspase inhibitors (blue line). 
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Table 1: List of proteins in Figured] and their descriptions. 



symbol 


description 


p53n 


nuclear p53 


p53T 


stress-induced p53 for apoptotic activation 


p53C 


p53 in cytoplasm 


p53CU 


ubiquitinated p53 in cytoplasm 


p53NU 


ubiquitinated p53 in nucleus 


p53U 


polyubiquitinated p53 for proteasomal decay 


p53m 


mitochondrial p53 


p53no 


p53n oligomer for transcription activation of Mdm2 


p53To 


p53T oligomer for transcription activation pro-apoptotic PUMA/Bax 


Mdm2-mRNA 


mRNA of Mdm2 


Mdm2n 


Mdm2/E3 ligase proteins in nucleus 


Mdm2C 


Mdm2/E3 ligase proteins in cytoplasm 


PUMA 


p53-upregulated modulator of apoptosis 


Bcl-2 


anti-apoptotic Bcl-2 family proteins 


Bax 


proapoptotic Bcl-2 family members 


Bax„ 


oligomerized pro-apoptotic Bcl-2 family proteins 


Baxm 


proapoptotic Bcl-2 members on the mitochondria 


Bid 


BH3-only pro-apoptotic proteins 


tBid 


truncated Bid 


cytc 


cytoplasmic cytochrome c 


Smac 


cytoplasmic Smac protein 


Apaf 


Apaf-1 pro-apoptotic protein 


XIAP 


IAP-family proteins 


apop 


apoptosome 


C9 


caspase 9 


C3 


procaspase 3 


C3* 


active caspase 3 


C3U 


ubiquitinated caspase 3^ 


C8 


caspase 8 



Table 2: Radiation dose dependent kinetic parameters for the p53 module. Ku and « 12 
are also dose dependent, but their values do not affect the observed, dose-dependent 
p53 oscillatory behavior and are noted in the text and the figures. 





No radiation (N-IR) 


Low radiation (L-IR) 


High radiation (H-IR) 




fCT 4 


fCT 4 


fcr 3 


K 5 (S" 1 ) 


2 x fCT 4 


2 x fCT 4 


5 x fcr 3 




fCT 4 


fCT 4 


5 x fCT 4 


k+ (/iM/s) 





2 x fCT 7 


5 x fcr 7 


^20 (S _1 ) 


fcr 4 


fCT 4 


fCT 5 


«23 (S _1 ) 


fCT 4 


fCT 4 


fCT 5 
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Table 3: Ki netic parameters of the p53 module that are unaffected b y radiation. 
Formation reactions: /iMs -1 

K = irr 6 

Unimolecular reactions: s~ 4 



k 3 


= 10~ 4 


k 4 = 


10- 1 


k 6 = 5x 10" 4 


K 


r = o 


K 8 = 


10- 3 


k 9 =0 


Kg 


= 10- 3 


«10 = 


10" 4 


K U = 10" 5 


^13 


= 10- 4 


«14 = 


= 0.1 


fc- = 5 x 10- 5 


«16 


= 10~ 2 


«17 = 


= 0.1 


Kis = io~ 2 


«19 


= 10" 4 


«21 = 


10- 3 


K22 = 10~ 3 


«24 


= icr 3 









Bimolecular reactions: //M s 1 
k$ = 0.1 kg = 0.1 

Special (cooperative) reactions 
k+ = 10 /iM^s- 1 fc+ = 10 //MV 1 
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T able 4: Kinetic parameters downstream of the p53 modu le. 
Formation reactions: /zMs -1 
1(T 5 for Bax, Bcl-2, Bid 
10~ 4 for Apaf, C9, C3, XIAP, Smacm, cytc 
Unimolecular reactions: units s _1 

K 25 = 10- 3 K 26 = 10- 3 k 27 = IO- 3 

fc" = IO- 3 fc+ = 5 X ID" 3 k 29 = 5 X ID" 3 

fc7 n = 10-3 fc+ = 5 x 10-3 L", = 5 x 10-3 



^36 = 


ID" 3 


k 37 = 


10-3 


k 40 — 


10- 


k 41 = 


10- 5 


k 42 = 


10~ 4 


k 43 = 


10- 


K43 


= 1 


k 44 = 


ID" 3 


K 44 


= 1 




10-3 


K 45 - 


= 1 


k 46 = 


10- 


k 47 = 


10-3 


K 47 = 


= 0.1 


k 4S = 


10- 


K^s = 


= 0.1 


k 49 = 


ID" 3 







10-3; PUMA, tBid, Bid, C3, C9, XIAP, Apaf 
lO- 3 : cytcm, Smacm, Smac, C3*, C3U, C9U 
lO- 3 : Bax, Bax*, Bcl2 
5 x 10- 4 : Mdm2, Mdm2C 
lO" 2 : cytC 



Bimolecular reactions: 


: /iM-^- 1 




k+ = 3 x lO" 3 


fc+ = 0.1 


fc+ = 0.6 


fc+ = 1.0 


k 31 = 2 x lO" 2 


/t 33 = 0.5 


K34 = 5 x lO -2 


k 35 = 0.5 


He = 0-6 


kf 7 = 0.1 


K 38 = 10 


/t 39 = 10 


k 4 Q = 0.3 


4 = 3 x lO" 2 


kf 3 = 3x lO" 2 


^44 = 0.3 


kt 5 = 0.3 


fc+ = 3 x lO" 2 


&47 = 1 


k 48 — 1 


^ = 0.1 



Special (cooperative) reactions 
fc+ = 10 5 //M^s" 1 
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